Attosecond tracking of light absorption and refraction in fullerenes 
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The collective response of matter is ubiquitous and widely exploited, e.g. in plasmonic, optical and electronic 
devices. Here we trace on an attosecond time scale the birth of collective excitations in a finite system and find 
distinct new features in this regime. Combining quantum chemical computation with quantum kinetic methods 
we calculate the time-dependent light absorption and refraction in fullerene that serve as indicators for the 
emergence of collective modes. We explain the numerically calculated novel transient features by an analytical 
model and point out the relevance for ultra-fast photonic and electronic applications. A scheme is proposed to 
measure the predicted effects via the emergent attosecond metrology. 
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I. INTRODUCTION 

The last decade has witnessed the emergence of the at- 
tosecond science opening a window on sub-femtosecond pro- 
cesses that take place in atoms and molecules (see 
Refs. [5-8] for reviews). Increased interest is currently fo- 
cused on many-body effects and condensed matter systems 
ll9i fTD where many degrees of freedom may interfere. The 
hallmark of extended systems is the collective, dielectric lin- 
ear response [ 12 1 that determines for example the light propa- 
gation lfl~3l n~4l and the energy and momentum loss of travers- 
ing charged particles lfl5"l[T6l . On a fundamental level, the di- 
electric response describes how the particles cooperatively act 
as to screen the interparticle Coulomb interaction. Obviously, 
this collective motion builds up on a time scale on which the 
effective interparticle interaction changes qualitatively (even 
in sign) [ 17 1. This was impressively demonstrated by terahertz 
spectroscopy for a semiconductor-based electron-hole plasma 
|[T8l evidencing that the electron-electron interaction develops 
from its unscreened to a fully screened form on a time scale 
of the order of the inverse plasma frequency (several 10~ 14 
seconds). For a finite system the situation is qualitatively dif- 
ferent: Here quantum size and topology effects bring about 
new features. E.g., as detailed below a spherical shell system 
exhibits two coupled plasmon modes that, while interfering 
transiently, evolve to two separate features in the long-time 
limit. In addition, these modes lie energetically in the ultravi- 
olet (UV) or extreme ultraviolet (XUV) energy range. Hence, 
tracing their evolution requires sub-femtosecond resolution, 
which is within reach experimentally |[]}|8). 

To be definite we concentrate here on the attosecond (as) 
evolution of the dielectric response in finite systems as moni- 
tored by the light absorption and refraction in fullerenes. The 
time-non-resolved (we call it hereafter steady-state) absorp- 
tion of the fullerene lfl9H2Tll . in our case Cqq, can be ac- 
cessed by analyzing the absorbed components from an im- 
pinging moderate intensity, broadband linear-polarized pulse. 
To monitor the attosecond buildup in absorption, we suggest 



* andrey.moskalenko@physik.uni-halle.de Also at A.F. Ioffe Physico- 
Technical Institute, 194021 St. Petersburg, Russia 



that at a time t = an XUV attosecond pulse ionizes Cgo 
and the change in the absorption at a time delay r D thereafter 
is recorded, as sketched in Fig. [T] For negative r D we obtain 
the absorption of Cgo and for large (few femtoseconds) r D we 
expect the steady-state absorption of C^q. Equilibration in the 
limit t oo proceeds via electronic, plasmonic and finally 
ionic channels, each having a distinct time scale; the ions are 
frozen at the as timescale. We focus on the plasma frequency 
(w p ) regime where the absorption and refraction show strong 
modulations l22l . The change of cj p for t > is evident 
from a qualitative consideration of the square root dependence 
on the density: The valence band of C$q accommodates 240 
electrons (originating from the 2s and 2p states of each car- 
bon atom) that form 180 er-type and 60 7r-type molecular or- 
bitals 1 20]. The highest occupied molecular orbital (HOMO) 
is five-fold degenerate. Its complete depletion by the XUV 
pulse reduces the number of particles by 4%, or red-shifts the 
collective response by 2%. As theoretically quantified below, 
even a single-photoionization event has a significant impact 
on the plasmon dynamics. 



II. THEORETICAL FORMULATION 

Ceo has a size below 1 nm [21], its electronic and opti- 
cal properties for t < and t — > oo are well-documented 
||20ll2TI . Also Cgo is readily available, very stable, and has an 
ionization potential of «7.5 eV and an electron affinity of 2 
eV It also exists in a well-studied crystalline form (fullerite) 
ll23ll . The collective response is marked by a giant plasmon 
resonance at s» 22 eV ll24ll that was investigated experimen- 
tally in the solid 11251 l26l and the gas phase [27]. Consider- 
able efforts were devoted to clarify quantitatively the experi- 
ments [28, 29 1 . It was not until recently however that the ex- 
istence of a second plasmon peak at higher energy (sa 39 eV) 
but with a lower oscillator strength was experimentally con- 
firmed l30ll3TI . Theoretically, the presence of two collective 
modes follows from a classical dielectric shell model P2~U36"1 . 
The quantitative description, as we are seeking here, is a chal- 
lenging task, however. Quantum mechanical approaches uti- 
lize either the tight-binding (TB) model for the valence elec- 
trons |24) or the jellium shell model lf3Tl [37T - i40l . The linear 
response to an external electric field is calculated either using 
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FIG. 1. (Color online) A schematic of the considered setup. An 
attosecond XUV pulse ionizes a fullerene. After a time delay r D the 
time-dependent absorption of the sample is studied by analyzing the 
absorption of broadband, linear-polarized pulse. 



the random phase approximation (RPA) for the polarization 
propagator (x) (in the case of the TB [24J and some of the 
density-functional-based (DFT) models |38 , 39 1), or by im- 
plementing the time-dependent DFT (TDDFT) OTl [37l l40l . 
Until now TDDFT seems to deliver the best agreement with 
experiments, however, the theory describes correctly either 
the position of the high energy plasmon resonance or, with 
the help of an adjustable parameter, the position of the lower 
energy plasmon resonance IT3T1 l40l . Here we develop an 
approach that rests on three steps: 1. We utilize quantum 
chemical ab-initio techniques to capture accurately the sta- 
tionary, single-particle electronic structure. 2. These states 
are expressed in a basis appropriate for many-body calcula- 
tions which we perform within the linear response theory (i.e., 
the random phase approximation) to obtain the steady-state 
collective response. 3. With both quantities at hand we per- 
form quantum kinetic calculations within the density matrix 
formalism to obtain the time and the frequency dependent po- 
larizability a(E, t = r n ). The imaginary and the real part of 
a(E, t) deliver then respectively the time-dependent absorp- 
tion and refraction properties of the sample |22|. 

We showed recently II4T1443I that the wave function <& e (r) 
of the valence band electrons with an energy e = e n i, as 
obtained from first principle calculations, is expressible to a 
good approximation as a product of a radial part R n (r) (hav- 
ing n — 1 nodes) and an angular part described by spheri- 
cal harmonics Yj m (f2) with l,m being the orbital and mag- 
netic quantum numbers (r = (r, tt) identifies the electron 
position, cf. Ref. 04)) l38l ATI . This procedure is shown 
ED to be valid for C 48 N 12 , B 80 , C 60 , C 240 , C 540 , C 20 H 20 , 
and Au7 2 . Thus, the theory presented below is readily ap- 
plicable to these systems; for the sake of clarity the discus- 
sion is restricted to Ceo, however. For Ceo the two occupied 
radial subbands £ n= \ i and £ n=2 ; are separated by approxi- 
mately 17.5 eV. The HOMO-LUMO gap as determined by our 
ab-initio calculations (E g « 5.5 eV) as well as some struc- 
tural information are encapsulated in the energy spectrum e n i 
(cf. fl4l ). These data allow us to perform the mapping of 
the unperturbed system to the single-particle Hamiltonian Hq 
with the states (f> a = <p n i m = R n (r)Yi m {VL) and the spectrum 
e n i- These single particle states are then taken as the basis to 
express the density operator p. Our theoretical description of 



the charge dynamics is based on the resulting density matrix 
P- 

We are interested in two types of external perturbations 
r/ext mat trigger the evolution: (i) a broad band, low-intensity 
light pulse probing the plasmonic response and (ii) a sudden 
change in the population upon a single ionization of Cgo by 
the attosecond XUV pulse. Formally, the latter change is de- 
scribed by a time-dependent occupation function f n i(t). To 
simulate the plasmonic response, we employ the Heisenberg 
equation of motion for the density matrix p, the mean-field 
approximation and the linear response to the external driv- 
ing. The relaxation within a time r due to collisions we treat 
within the particle-conserving relaxation time approximation 
P31l . This approach has been successfully tested for a variety 
of systems [18, 46-49 1. It should be emphasized, however, 
that we are interested in the dynamics at times much shorter 
than r. The equation of motion for p reads as 



dp i 
di + h 



H + V ind + V c *\p 



-I.e. 



(1) 



The induced potential V lnd has to be determined self- 
consistently from the induced charge density as derived from 
the change in the density operator p [to find p we need V" md 
in Eq. ([TJ]. p l c - is the local equilibrium density operator. The 
corresponding locally relaxed density p 1 °- (r, t) is the distri- 
bution that, at any given instant, would be in equilibrium in 
the presence of V ext and V while satisfying the charge 
conservation II45L Technically, we express p in the basis 
{4>a(r)} an d expand to first order around the equilibrium, i.e. 
p a p(t) = (a\p(t)\(3) = f a (t)S aP + 6 PaP (t); p x a %t) = 
(a\p le -(t)\f3) = f a (t)5 afi + n° a0 (t)5 Pa0 (t), where n° a0 (t) = 

an< i i s the matrix corresponding to the local 
chemical potential (full details are given in Appendix|A]). 

Application of an external linearly-polarized electric field 
with an amplitude £ , polarization direction e, and fre- 
quency uj, that corresponds to V cxt (t) = — er ■ e£oe~ lult 
(e is the electron charge), leads to the change of the den- 
sity matrix 8 pit) which in linear response is governed by 
the two-times response function x(tjt'), i.e. 8p(t) = 
J_ dt' £')V cxt (i') . The induced dipole moment derives 

from P(t) = Tr [er8p{t)\. For spherical molecules, reduc- 
ing all quantities to their radial components, we find (see Ap- 
pendix |E| 



4tt 



e 2 £ / dt' r T xit,t')re 



P(t) 



where the elements of r are r nn i = r n > n — (n\r\n'). Intro- 
ducing the dimensionless, Fourier-transformed response func- 
tion (f = r/rn, r is the fullerene average radius and eo is the 
vacuum permittivity) 



z(ui,t) 



/ dre^x(t,t-T)r 
Jo 



(2) 



we find P(w; t) = -4-Kr^e Q £ e~ lut r T z(w, t) . The time and 
frequency-dependent dipolar polarizability is then calculated 
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FIG. 2. (Color online) Energy dependence of the (a) real and (b) imaginary parts of the dipolar polarizability a(E) of Ceo in the long-time limit. 
Insets to (a) show the calculated spatial densities of the plasmonic modes, corresponding to the peaks in the spectrum, in any plane crossing the 
center of the molecule and being parallel to the polarization of the external electric field. Dotted lines result from the approximation including 
only intraband collective excitations. Inset to (b) shows the difference in the spectra between Ceo and C^ for the case of refraction. 



numerically as 



a(u, t) 



3^T 



-rXr 



z(u),t) 



(3) 



(in SI units, asi = ifteoa). The explicit equation for the 
evolution of z(u, t), and therefore of a(uj, t), following from 
Eq. ([Tji we included in Appendix [A] 

For an insight into the numerical results we construct an 
analytical model based on the following: For Cqo (and gen- 
erally for „spherical" molecules) we find that the intraband 
components z n (w,t) and 222(^5 1) of the response function 
are hardly influenced by the interband components 212(10, t) 
and Z2i(u>, t) (details are in Appendix [d]). In this approxima- 
tion, considering only the intraband terms we infer a system 
of two coupled equations for Zh(uj, t) and -222(^7 i) 



1 \ 2 2 W 

2iuj \ Of + u> F _ — u —1 — 

T J ' T 

n' = l,2 



(4) 



where the stationary matrix g is explicitly defined in Ap- 
pendix IB] and Wp „ = 1/h 2 x e 2 \f2N n / (§ire a r ) x fiwF,n- 
Here iv^is the number of electrons in the n-th radial sub- 
band. For Cgo in the equilibrium we have N\ = 180 and 
N2 = 60. hjj-p,n are the energy distances between the highest 
occupied state in the ?i-th radial subband and the unoccupied 
state in the same subband having the value of I that is greater 
by one. On the same approximation level Eq. ([3) simplifies 
to a(u>, t) = — Tq [z%i(u), t) + Z22(w, f)] meaning that the in- 
terband components can be then neglected when calculating 
a(u!,t). In such a consideration the dynamics of the electron 
density in two different radial channels are still coupled via 
the induced fields like in concentric nanoshells 1501. 



III. STEADY-STATE COLLECTIVE RESPONSE 

Absorption and refraction spectra of the Cgo and Cg~ in the 
long-time limit (i.e. r D < 0) are shown in Fig. [2] The re- 



sponses of C^g and Ceo are quite similar [cf. inset in Fig.[2|b)] 
which is in agreement with the existing theory and experiment 
OT1 l40l . Our Ceo calculations exhibit two peaks at 23.2 eV 
and at 39.6 eV in the absorption spectrum which correspond 
to plasmon excitations and agree very well with experiment 
||3T| . The width of the measured main plasmon resonance is 
reported to be « 5 eV which sets the scale for the relaxation 
time r. With a(E) = ^Ima(£) E2, where c is the speed 
of light, we find the corresponding absorption cross-sections 
a(E) to be 4.1 x 10 3 Mb for the main peak and 9.4 x 10 2 Mb 
for the higher energy peak. This yields relative heights of the 
peaks that are in qualitative agreement with the experiments 
and with the time-dependent DFT results ll37l l40l . There is 
also a weak low energy contribution at m 9 eV. Such struc- 
tures in the absorption have been attributed to single-electron 
transitions [20 28 1. Our calculated resonance has a collective 
character and is energetically in the proximity of the single- 
electron excitations making it difficult for an observation via 
absorption. As well-established [22 1, the peaks in the absorp- 
tion (refraction) are symmetric (antisymmetric) with respect 
to reflection at the respective central frequency. 

The approximation including only intraband collective ex- 
citations, i.e. considering the stationary solution of Eq. Q, 
gives a very good description of the main peak in the spec- 
trum (see Fig.[2]i, but misses the higher energy plasmon. The 
last is caused by interband collective excitations (i.e. gov- 
erned by Z12 and 221) that are still influenced strongly by the 
intraband response. The interband excitations shift slightly 
the main plasmon resonance to lower energies and enhance 
its strength. The /-sum rule ||39l IBTl l52l for the absorption is 
fulfilled only approximately, particularly because higher en- 
ergy states in the continuum and close to it are not included 
in our model. However, these states contribute mostly to the 
background of the spectrum. For a further insight we show the 
calculated spatial distributions of the induced charge densities 
at the maximal absorption [insets in Fig.[2ja)]. The oscilla- 
tion amplitude depends linearly on the external field strength. 
The interband character of the the higher energy plasmon ex- 
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FIG. 3. (Color online) (a) Difference between the imaginary part of the dipolar polarizability lma(E,t) after removing one electron from 
the highest occupied state of Ceo in the second (upper) branch at t = and its long-time limit for Cg , i.e. Ima c + (E). (b) The same 

dependence for the real part change of a(E, t). (c) Time dependence of the imaginary part of the dipolar polarizability Im a(E, t) of the 
excited fullerene molecule at fixed energy values, E — 24.6 eV and E — 38.8 eV, taken close to the main and higher energy resonance 
positions of the stationary response, respectively. Time intervals T = nh/E indicated by the horizontal double-arrow lines are close to the 
observed oscillation periods for each of the dependencies: T = 84 as and T = 53 as , respectively. 



citation is evident. Its radial structure is determined by the 
product Ri (r)i?2 (r) and has therefore one node located in the 
middle of the fullerene cage. The angular dependence of the 
modes in and out of the depicted plane is governed by their 
dipolar character. The radial structure of the density oscilla- 
tion at the main plasmon peak shows a constructive superposi- 
tion of the radial electronic radial density of the first |i?i(r)| 2 
and the second |i?2(r)| 2 radial channels that oscillate in phase. 
In contrast, we observe an out-of -phase oscillation of the den- 
sities in different channels at 9 eV. Both cases are collective 
intraband electronic excitations. 



IV. TRANSIENT DYNAMICS 

Having assessed successfully the steady-state response we 
move on to the attosecond transient dynamics. We focus on 
the setup shown in Fig.[T] At t = the XUV pulse with en- 
ergy uj oj p and a duration < 100 as changes the population 
swiftly in Cgo by ejecting one electron. To access the transient 
absorption or refraction let us assume the photoelectron is 
emitted at t = from the second radial band leaving Cg (we 
find similar results for multiple XUV ionization). For all time 
moments before the ionization the time-dependent response 
a(E,t < 0) is equal to its steady-state value ac 60 (E) for 
C6o- At t > we solve for the dynamics of a(E, t), which for 
long times approaches the steady-state value a c + (E) for C^q. 
The results in Figs.[3ja) and[3|b) illustrate the change in the 
imaginary (absorption) and real (refraction) parts of the po- 
larizability by evaluating Aa(E, t) = ac 60 (E, t) — a c + (E). 
The transient dynamics shows a rich structure evolving over 
approximately 2 fs. For 1 fs the remainder of the vanish- 
ing difference between a(E. t) and a n + (E) is still visible 

_ _ Wo 

in Figs.pla) and Bib). For t < 100 as the response is basi- 
cally determined by that of Ceo, i.e. it takes around 100 as 
for Cg" to start responding collectively and it attains its full, 
steady state response for t > 1 fs. The reason of this iner- 
tia is the finite mass of the carriers. Three marked general 
transient features can be distilled from Fig. [3] The dispersive 



hump starting around 200 as is due to the sudden change in 
the population (and hence the wide frequency perturbation) 
brought about by the shortness of the XUV pulse. Further- 
more, in addition to the relaxation to the steady-state response 
we observe an oscillatory behavior in time with a certain fre- 
quency. This behavior is most obvious around cj p (cf. Pig . [3j . 
The origin of these oscillations is comprehensible from a con- 
sideration of the response of the fullerene molecule in the 
lowest approximation: As we demonstrated by Eq. Q the re- 
sponse resembles two coupled damped harmonic oscillators. 
For photon frequencies u close to the main peak position the 
model can even be further simplified and the spectra are de- 
termined approximately by the dynamics of a single driven 
damped harmonic oscillator with the frequency correspond- 
ing to the main peak position. Hence, the feature in the re- 
sponse we can now understand from the known properties of 
the driven, damped harmonic motion. On a short time scale 
this dynamics is governed by a combined decay and oscilla- 
tions with a frequency around 2cj, which is clearly observed in 
the full-fledge calculations shown in Fig. [3jc) for Im a(E, t) 
at E = 24.6 eV (time period corresponding to 2cj is indicated 
by a horizontal double-arrow line for comparison). The dy- 
namics of lma(E,t) for fixed energy E = 38.8 eV in the 
vicinity of the higher energy peak is contributed to by more 
frequencies because the higher energy plasmon is strongly in- 
fluenced by the main plasmon. However, oscillations corre- 
sponding to 2lj are also seen. 



V. SUMMARY 

We developed a framework for the attosecond collective re- 
sponse in finite systems with spherical symmetry and applied 
it to fullerenes. The predicted marked features in the transient 
absorption and refraction should be accessible with current at- 
tosecond metrology. The discovered attosecond dynamics in 
the optical response points to new opportunities for optoelec- 
tronic devices at the sub-femtosecond time scale. 
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Here e n i m are the single particle energies (cf. Ref. B4| ). The 
matrix elements of the Coulomb potential are cast as 



V nlm , M = (nlm\V lnd + V^\n'l'm') 



ind i y/cxt 

nlm.n'l'm 1 3 



V, n, ,+V. 
v nlm.n' I'm' 1 v •* 



(A7) 



Appendix A: Derivation of the linear response equations 

The transient time and energy-dependent polarizability 
a(t, E) is governed, to a first order (linear response), by the 
two-time response function if). For the determination of 
x(t, t') we setup a calculational scheme based on the solution 
of the Heisenberg equation of motion for the density matrix p 
under the influence of the external field V rcxt in the mean-field 
approximation that leads to Eq. ([T|. In the basis {4> a } we find 
the density matrix elements as p n im,n'l'm' — (nlm\p\n'l'rn'). 
From the change in the density matrix 5p n im,n'i'm' we calcu- 
late the change in the charge density 



where 



and 



T^CXt 



nlmm'l'm' 



E dZw'Vffir' (A8) 



t ?-ext,n7 
V l"m" 



r 2 dr Snn> (r)V^„(r) 



(A9) 



Upon solving the Poisson equation (see Appendix |B]> we find 
for the matrix elements of the potential V md (r) = e$(r) in- 
duced by the change in the charge density 8n(r) the following 
expression 



Sn(r) = ^2Sni m (r)Yi 



ml", <+>! , 



in i 



where 



'(r)Sn 



Im ' 



s n 'n"(r) = R n ,{r)R n n{r) , 



(Al) 



(A2) 



(A3) 



An™ " 



htl'm'V'rn'' 



2_j y\m,l"m"^Pnlm,n'l'm' > (A4) 
I'm' ,l"m" 



'(2/' + l)(2J" + 
4tt(2Z + 1) 



CIO /~ilm 
V 0,1" 0° I'm' ,1" 



(A5) 



Cfc; i„ m /, are Clebsch-Gordan coefficients. The change in 

the charge density determines the induced potential V md lead- 
ing thus to a self-consistent procedure. 

In the basis {(f) a } the Heisenberg equation of motion for the 
density matrix ([TJ reads: 



&Pnlm,n'l' 



i 

h 

i y 

h 4? 

n"l"m' 



Pnlm,n f l'm' 



{,Pn"l" m" ,n'l'm' ^ / n/m,n"Z // m" 
Pnlm,n"l" m" ^n" I" m" .n'Z'm' ) 



r)-- c - — n 

Pnlm.n'l'm' Pnlm.n'l'in' 



(A6) 



ind = \^ \^ lm I" At)™ 1 '™ 2 

nlmm'l'rn' / , / , vl' m' J" m" ynn'ri\ri2 I" m" 

en — — 

rii«2 l"m" 



A10) 



V. 



The coefficients g nn , n n2 are defined in Appendix B by 
Eq. (|BT3). 

We are interested in the linear response and expand there- 
fore the density matrix around the equilibrium distribution 



J a 



Pnlin,n'V m' ~ fnlm ^ Prdm ,n' V 



(All) 



Thus, 6pnlm,n'l'm' is a small perturbation. The elements of 
the local equilibrium density matrix including the lowest order 
correction with respect to the equilibrium density matrix are 
given then by B31 



Pnlm,n'l'm' ~ fnlm^nn' Sa'S mn 

fO _ fO 
J nlm J n'l'r 



& Pnlm.n'l'r 



(A12) 



Here the components of the chemical potential can be ex- 
pressed as 



l f 7n , .n"l"m' 



= Vi/!!, m '„?-W'"" • (AO) 

/ i a l"m" ,lm ^lm v ' 



We require density conservation, meaning that n r {, 
n)'^ ,n n which in turn leads to 



"lm 



where 



lm,lm 



Sn&" = Y M, n 'f- 5tf>" , (A14) 



fO _ fO 

EJm< ir rn ' n'l'm' Jn"l"m" 

yim,l"Tn" yi"m" Jrh £ , 



I'ra'V'm" 



£n"l" m" 

(A15) 

In general, SpPJ 1 is determined from Eq. ( |A14| i by solv- 
ing a system of linear equations for each pair of indexes 
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n',n". This step of the calculation is significantly simpli- where 
fied by neglecting the energy splitting of the multiplet with 
the same I, i.e. by writing £ n i m — £ n i m ' = £ni an d there- 
fore f° lm = f° lm , = /",, which is a good approximation for 
"spherical" molecules such as Ceo, as confirmed by the ab- 
initio calculations. Adopting this approximation we conclude and 
that (cf. Appendix [C| 



Mr 



VI" £n ' 1 ' 



fO _ -fO 
Jn'V Jn"l" 



£n'V £n"l'' 



l (2/' + l)(2J" + l) 



4tt(2Z 



\yvo,l"o\ 



(A17) 



(A18) 



Gathering the informations in Eq. ( | A 1 3] > and Eq. ( |A16[ ) we 
can state 



I'm' ,n"l"m' 



= Yytr', r 

/ j a l"m" 



1 



— 6n?. 



Irh ]\/[ n ' n " I 



(A19) 



5(i 



n n 
I in 



M n'n 



— <KtT 



rw 

at 



The solution of this equation is found as 



(A16) 



From Eq. \A6) we deduce then the following determining 
equation for the change in the density matrix 



i(£ n i - e n 'l>) + 



5 Pnlm,n'l f 



1 Jnl Jn'l' r„ 

° Pnlm ,n' V 

T E n l — E n 'l> 



5pnlm,n>l>m'(t) = % - j dt 1 e^-^H] <*"*'> | [/°,(*0 - VW.n'IW (0 

~r t u f^nlm.n' V m \ L J 



(A20) 



(A21) 



Inserting Eqs. ( |A10| >,( |A8] > and ( | A 1 9| > into Eq. ( |A21[ ) and summing the left and right hand sides weighted by y\™i" m n as in 
Eq. we infer that 



An™ " 
0n lm 



At' nf'«" \t,H)v£' n ' w V) + / &t> n^"{t,t') e r J2 9 l n > n », ni nJnZ n2 V) 



+ dt'ir 1 (t,t')5n? m n (f), 



(A22) 



where 



nr'""(^*') = |E^^" e[l(e "' ! '" e "" , " )+ ^ ](t '" 4) [fn'At 1 ) - M)] 



(A23) 



i'i" 



if n "(t,t') = l — Y KLj^'i'-^)H]it'-t) fi'i>(f) (A24) 

tM 1 ™ ^ £ ™'*' _ £ ""i" 

This is the central integral equation, which has to be solved to arrive at the time-dependent density change induced by a perturbing 
external electric field. 

We introduce the response function x l nn ' ni n 2 (*> O as 



Snf m n "(t)= f dt> J2x l n'n»,n in ^W™' nin2 (t') 
J — oo „ 



(A25) 



Here we consider the perturbation by the external electric field in the dipole approximation. Therefore the calculation can be 
limited to the case Z = 0. In order to simplify notations, in the following we omit the corresponding index by all variables and 
coefficients. 
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Inserting the definition ( |A25[> into Eq. ( |A22| i, changing the order of integration in the second term on the right hand side, 
comparing again with Eq. ( |A25| > and making use of the dipole approximation, we get the following integral equation for the 
response function: 

n 1 n 2 

r* P 2 , n ^ 
_i_ / Af" —U n n (t t") n , ,, v (t" t')r 

^ I ut/ A± {"i I / t yn'n / \nin 2 Anin 2 , 713714 X 1 * i L )' 713714 (A26) 

^ nin2,n3ri4 

4- f At" T n ' n " (t t") V y / « If" t')r 

*'*' n t n 2 

where the matrix elements f nn > are given by rv„» = i Jr 2 dr s n ' n "(r)r and are calculated by us using the radial functions 
for Ceo from our ab-initio calculations. Here we have also taken into account that in the dipole approximation the identity 
V cxt ' nin2 f n r n ir = V cxt ' n n r ni n 2 holds. The corresponding time- and frequency-dependent response function is determined 
then by 



2 roo 

Z n ' n "(u},t) = / (iT e iuT ^ Xn'n",nin 2 (t,t-T)r nirl2 , (A27) 

£ o Jo 



where g = l/(3ro) and eo is the vacuum permittivity. Here tq 6.745a B is the average radius of the C$o atomic cage, where 
a B is the Bohr radius. After inserting Eq. ( |A26[ > into Eq. ( |A27[ ) we arrive at the following system of integral equations for each 
frequency value w: 



J — OO 

+ f dt'e^^W^^'^t') 9n'n", ni n 2 Z ni n 2 ^,t') (A28) 

n u n 2 

+ ( dt'e^-^I^" (t, t')z n , n ,{u, t') , n', n" = 1, 2 ; 



Where g n 'n" ,nin 2 — QWn" ,nin 2 / ' 9 an d 

2 

W n ' n "(t,t') = ( -9-n n ' n "(t 1 t') 



(A29) 



Her e U n ' n " and I n ' n " ' (t,f) are given by Eq. ( |A23] > and 
Eq. ( A24[ >, respectively, with 1 = 1. The system of equations 
( |A28[ > is solved numerically. 

It can be shown (see section[E]) that the time- and frequency- 
dependent dipolar polarizability can be written as 



ol{uj, t) 



^ ^ ^n' ' n" %n' n" S ^) 



(A30) 



Appendix B: Calculation of the induced potential and its matrix 
elements 

Here we determine the electrostatic potential <£>(r) induced 
by the change of the density in the spherical layer n(r). The 
Poisson equation 



A$(r) 



-n(r) 



(Bl) 



has then to be solved (SI units are used in this work). Its 
solution can be written using the Green's function G(r,f / ) 



$(r) = -— / G(f,r )n(r)A 3 r 



where 



1 



47r I f — r* I 



(B2) 



(B3) 



Using the spherical geometry of the problem we apply the 
following decomposition for the Green's function: 



OO I 



1 A * 



i=o m=—l 



(B4) 

where r> = max{r, r'} and r< = min{r, r'}. Then we get 
the following expression for the potential 



$(r) 



OO I 



E 



1 rl 



1=0 m=-l 



21 + 1 ri+i 



1", 



> 

xYr m (0',«£')n(r<) 



(B5) 
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To proceed further we decompose n{f ) in spherical har- 
monics 



lir, 9,(j)) = y^ j ni m (r)Yi m (9, <j)), 



(B6) 



where the radial-dependent angular components of the density 
can be found as 

ni m (r) = ^2 s n>n"(r)rip r £ , (B7) 

n'n" 



with 



s n 'n"(r) = R n '(r)R n "{r) , 



(B8) 



and finally expressed via the components of the density matrix 
Pnim,n'Vm' because 



E 



Viral" m" Pnlm,n'l'm' 



(B9) 



I'm' M' 



After inserting Eq. < |B6| > with Eq. < |B7| > into Eq. < |B5| > we get 



$ ( r >-EEE rC'ir'(r)>WM) , (BIO) 
£ o , — , 

nn' I — U m— — l 



where 



Il m (r) 



1 



21 + 1 



dr' r' 2 ^s 



r>) 



(Bll) 



The determination of the induced electrostatic potential ( |B10| > 
entails the knowledge of the elements of the density matrix 

We calculated the matrix elements 

(n'l'm'\<f>(r)\n"l"m") = $> n 'l'm',n"i"m" of the induced 
potential as 



n'l'm' ,n"l"m" 



£q y l"m",lm yn n > n i n 2 Ifh. ' 



nm 2 iA 



where y 1 ™ n is given by Eq. ( |A5| l and 



(B12) 



<?„v,w = /dr r'wWCM • (B13) 



The values of all matrix elements entering Eq. ( |B13| l are ob- 
tained numerically using the radial functions for Ceo from our 
ab-initio calculations. 

As explained above, for Ceo we may adopt the approxima- 
tion of a spherical layer with the width of the layer d being 
smaller than the sphere radius d « rj. This enables us to 
write 



9 n 'n" ,nina ~~ 9l (^nin 2 ^ 



n" !!„'„■ 



im ^9n'n" ,ni ni ) ' 

(B14) 



where 



9i = 



1 1 

27+1 ro' 



(B15) 



Sn'n'.mm = gn'n'.n.m ~ 1 and ^n'n",nm 2 < 1 so that 
the second term in Eq. ( |B14| i can be neglected to a the leading 

order of d/r^. 



Appendix C: Chemical potential calculation 

With e nlm = e nlm , = e ni and f nlm = f° lm , = jf°, 
Eq. ( |A15[ > can be written as 

fO _ fO 

T^fn'n" \ ^ Jn'V Jn"l" \ ^ i^' m ' 7/^ m ' 

lm,lra ~ 2-^1 e ,,, _ p „,„ 9lrn,l" ' m" it \n m ii Tjh ■ 
j, ,., c n'l' c n"L" . ,, 
v .1" m,m" 

(CI) 

The last sum on the right hand side of this equation we evalu- 
ate using the following property of the Clebsch-Gordan coef- 
ficients 



Es-tjm s-fj'm' 
iimi jimi jiroi 



j 2 m 2 ~ fijj'dmm 



(C2) 



(C3) 



Thereby we make use of the definition ( |A5[ ) and rewrite 

1l l ' m ' n l ' m ' = (l\m+fh l-m 1/ ~ m 

ylm,l"m" vyim" Irh > ' °V —tn' 'm" '»(' —m',l'' 

Summing over m' and m" and using Eq. (|C2j we find 

E = S ltJ 5 m ,*Kl n „ , (C4) 



where K vv , is given by Eq. ( |A18| >. With this finding Eq. ( |CT 
simplifies to 



Im.lm M ' < 



(C5) 



where M" n is defined by Eq. (jATTJi. Equation ( |A14| i can be 
then simplified to 

6nfc" = M? n "5^" , (C6) 
and after dividing this by M" n we get finally Eq. ( |A16[ ). 

Appendix D: Approximate stationary solutions 

Utilizing the spherical shell approximation (d/ro <C 1) we 
find g n n,mn2 ^ 1 anc l ^in 2 1 for ri\ ^ n 2 whereas 
9nn. mm ~ 1 an d ^mm ~ 1- As a consequence it follows 
from Eq. ( |A28[ ) that the intraband components zn(uj,t) and 
Z22 (w, £) of the response function are only weakly influenced 
by the interband components Zi2(u>,t) and Zi2(u),t). There- 
fore, we obtain an approximate solution by solving the equa- 
tion system ( |A28| > considering only the intraband terms. This 
leads to a system of two coupled equations for zu(lu, t) and 

The equilibrium population of levels is determined by the 
energetic order of the states. The expressions for W nn (t,t') 
and I nn (t, t') can thus be written approximately as 

W™{t, t')=^ n e^ Sln Mt> ~ t)] , (Dl) 



1 l,*> 



I nn (t, t') = -e - cos [u F , n (lf - t)} , 



(D2) 
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where 



ui 



e 2 N n 1 



WF,n - E g e m 67rrg 



(D3) 



iV„ is the number of electrons in the n-th radial subband. 
hiOF,n are me energy distances between the highest occupied 
state in the n-th radial subband and the unoccupied state in 
the same subband having the value of I that is greater by one. 
This energy distance can be written as 



where <? U)2 2 
D = 



= 522,11 and 

( 2 - i 2 2 

I W p, 2 522,22 T Wp,2 — U) 



T 



( W p,l5ll,l 



W p,l W p,2 (511,22)' 



w F 1 — w —1- 



(D10) 



The stationary frequency-dependent dipolar polarizability is 
given then by 



hid 



F,n 



1 2 



g ' 



(D4) 



where the first term on the right hand side determines the cor- 
responding energy gap between the highest occupied state and 
the lowest unoccupied state for N n free electrons in a thin 
spherical layer with the radius tq and E s is an additional band 
gap appearing for the real "spherical" molecule Cgo- 

With kernels determined by Eqs. ( |D1| | and ( |D2| i the sys- 
tem of integral equations for t) and z 22 (ui,t) can be 
reduced to two coupled ordinary differential equations, each 
of which describes the plasmon dynamics in the respective ra- 
dial band: 

d 2 , s A „. \ d , v 



+ ( w p,i5ii,ii+ w f,i- w2 -*^) zn(u,t) 

+ Wp )1 5ii 1 22^22(w,t) = -Wp j , 



(D5) 



At 2 



z 22 



(uj,t)+ f~-2iwj ^z 22 (u;,t) 



+ (^p,2522,22+wf. 2 -W 2 -i^ Z 22 (uj,t) 
+ 2522,11^11 = 



y p,2 



(D6) 



Here we made use of the fact that the corrections to fn = 
r 22 = 1 are of a second order in d/ro and can be neglected. 
Then the time- and frequency-dependent dipolar polarizability 
is expressed as 

a(u,t) = -rg [z n (uj,t) +z 22 (uj,t)] . (D7) 

The stationary solution of the system (RJ) can be found as 



Zii(w) 



1 

D 



W p,l W p,25ll,22 



Wp 2 522,22 



^F.2 



(D8) 



2:22(0;) = 



W p,l W p,25ll,22 



w p,i5ii,ii 



LU-p j — CO — 1- 



(D9) 



h(cj) = -rjj [zn(cj) + z 22 (w)] 



(Dll) 



The plasmon frequencies as retrieved from Eq. ( |D1 1[ ) de- 
rive from the zeros of denominator (for 1/r = 0). Since 

w fu - w f,2 < w p,i5n,n - Wp j2 522,22 we find two possi- 
ble frequencies w + and u;_ such that 



1 



where 



W P + ^2 ( W F,l w p,l5ll,ll +^F,2 w p,2522,22), (D12) 



2 , ,2 ~ 

P 

2 , ,2 ~ 



(2 2 ~ 2 2 - \ 

w F,l w p,2522,22 +^F,2 w p,lflll,llJ 



2 2 ~ 2 ~ 

^p = w p,l5ll,ll + W p2 522,22 



(D13) 



(D14) 



The amplitude of the peak centered around the frequency ui- 
has much smaller value because it contains terms of the or- 
der (w F1 - w| i2 )/(Wp !l 5ii,n - Wp i2 522,22)- From Eq. ( |D14| ) 
we calculate w + = 25.4 eV and u;_ = 8 eV. Comparing the 
results found so with those without making the above approxi- 
mations confirms the accuracy of the approximate expressions 
for the peak positions. 



Appendix E: Time- and frequency-dependent dipolar 
polarizability 

Given the time-dependent perturbation of the system den- 
sity Sn(r, t) the induced dipole moment P(t) is determined 
then by 



P(t) 



rSn(f, t) A 3 f 



(El) 



Inserting here Eqs. and \A2\ we find 



?{t) = e W( 



(E2) 

The components Sri^ (t) can be calculated once we have 
the components VfL' n n (t) of the external perturbation and 
the response function X l n ' n " mn 2 usm g Eq. ( |A25) >. 

If we consider the response of the system to an external 
spatially homogenous time-dependent electric field £(t), it is 
enough to consider only the dipolar components of the re- 
sponse function, i.e. the components with I = 0. We may 
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select the z-axis to be parallel to the electric field. Then the 
external potential reads: 

V c ^(f,t) = -ez£(t) . (E3) 

Considering the angular components of this potential we con- 
clude that only the component with / = 1 and m = does not 
vanish: 

V^(f,t) = -2^er£(t) . (E4) 

As a consequence in Eq. ( |E2[ > only terms with I = 1 and m = 
survive and we infer that the induced dipole P(t) moment 
is directed parallel to the electric field £(t). Its magnitude is 
given by 

4:7r 

P(t) = -—e 2 rl / dt f Vn"Xn'7 l " 1 T il n 2 (M / ) 

" J— OO / // 

n n ,niri2 

(E5) 

where 

f n 'n" = — (r 2 dr s n ' n »(r)r . (E6) 
For a perturbation with an electric field £(t) = £oe~' lu ' t at a 



particular frequency ui we derive using Eq. ( |A27) > the expres- 
sion 

P{uj;t) = -4Trr 3 Q e ^ TVn"*n'n"(MS>e _<W * ■ ( E7 > 

n'n" 

The time- and frequency-dependent dipolar polarizability 
asi(w, i) in SI units is then found to be equal to 

a S i(uJ,t) = -47re ro ^ r n > n » z n > n »(uj,t) . (E8) 

n'n" 

In literature, however, the polarizability is more frequently ex- 
pressed as polarizability volume (as it happens naturally when 
using Gauss units): a = j^tCksi- Using this definition we 
end up with 

a(w,t) = -Tq ^ r n >n"Z n 'n"(u,t) . (E9) 
n'n" 

We note that the factors 2 n 'n" t) are dimensionless. 

For d/ro -C 1 and considering only the lowest order con- 
tribution we conclude r n / n » w 6 n ' n" an d therefore 

a(w,t) w-rg^jB^w.i) . (E10) 

n 

With this expression we can describe the approximate re- 
sponse given only by the intraband excitations. 
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